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This paper is devoted to the study of the processes e + e~ — > J/^r) c , J/^rf c , ip'r/c, ip'T)' c within light 
cone formalism. It is shown that if one disregards the contribution of higher fock states, the twist- 
3 distribution amplitudes needed in the calculation can be unambiguously determined from the 
twist-2 distribution amplitudes and equations of motion. Using models of the twist-2 distribution 
amplitudes the cross sections of the processes under study have been calculated. The results of the 
calculation are in agreement with Belle and BaBar experiments. It is also shown that relativistic 
and radiative corrections to the cross sections play crucial role in the achievement of the agreement 
between the theory and experiments. The comparison of the results of this paper with the results 
obtained in other papers has been carried out. In particular, it is shown that the results of papers 
where relativistic and radiative corrections were calculated within NRQCD are overestimated by a 
>— ^ , factor of ~ 1.5. 
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PACS numbers: 12.38.-t, 12.38.Bx, 13. 66. Be, 

I. INTRODUCTION 



Double charmonium production at B-factorics is very interesting problem from theoretical point of view. At the 
beginning, the interest to this problem was caused by large discrepancy between theoretical predictions for the cross 

i 1 section of the process e + e~ — * J/vP^c [E 0, HI and it's first experimental measurement by Belle collaboration j4j. 

; Lately, large discrepancy between theory and experiment was found by Belle [j| and BaBar Q collaborations for the 
*p. processes e+e~ -> J/^r)' c , ip'r] Cl tj)'r)' c , J/^XcO, i>'XcO- 

O i' A number of attempts were made to explain this discrepancy. For instance, the authors of papers @, S 0] studied 
the possibility of admixture of e + e~ — > J I if) J I'll) events in the e + e~ — > J /^>rj c . Another attempt was to attribute the 
i ' discrepancy to lar ge r adiative corrections [1C|, [lH EH ■ One more very popular approach to the problem [ID, EH EE 
EH, E3) EM EM I20I |21| consisted in taking into account internal motion of the quark-antiquark pair in charmonium. 
t-H ■ Some other interesting points of view on this problem were proposed in papers [22|, [23|, [24| . 

Among many approaches to the resolution of this challenging problem one should mention nonrelativistic QCD 
(NRQCD) [25l | . Contrary to many other models of quarkonium production this approach allows one to improve 
systematically the accuracy of the calculation. Thus at the leading order approximation [j, H, HJ there is very large 
discrepancy between NRQCD predictions for the cross section of the process e + e _ — > J/^r/c and the experiments. 
However, after the inclusion of the one loop radiative corrections [ToL El| this discrepancy becomes smaller. At the 
next step, after taking into account radiative and relativistic corrections simultaneously [lEHlJj the problem can be 
resolved. From this example one can draw a conclusion, that relativistic and radiative corrections play very important 
role in the processes with charmonia production. Probably, in a similar way one can resolve many other puzzles with 
quarkonia production [28| . 

It should be noted that within NRQCD the discrepancy between the theory and experiments was resolved only 
for the process e + e~ — > J j^r\ c . At the same time the discrepancies were observed also for the processes e + e~ — > 
J/\f r ?7c, ip'ric, ip'r/'c, J/^XcOi ip'XcO- It is not clear is it possible to apply the same approach for these processes, since 
excited charmonia states are rather relativistic and the application of NRQCD to the production of such states is 
questionable. 

Another systematic approach to the study of hard exclusive processes is light cone formalism (LC) [2!|[3(|- Within 
this approach the amplitude of hard exclusive process can be separated into two parts. The first part is partons 
production at very small distances, which can be treated within perturbative QCD. The second part is the hadroniza- 
tion of the partons at larger distances. This part contains information about nonperturbative dynamic of the strong 
interactions. For hard exclusive processes it can be parameterized by process independent distribution amplitudes 
(DA), which can be considered as hadrons' wave functions at light like separation between the partons in the hadron. 
It should be noted that within LC one does not assume that the mesons are nonrelativistic. This approach can equally 
well be applied to the production of light and heavy mesons, if the DAs of the produced meson are known. For this 
reason, one can hope that within this approach one can study the production of excited charmonia states. 

The first attempts to describe the experimental results obtained at Belle and BaBar collaborations within LC were 
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done in papers [IH, 0, [3 H3] ■ There are two very important problems common for all these papers. The first one is 
that in these papers the renormalization group evolution of the DAs was disregarded. What is very important since 
the evolution of the DAs takes into account very important part of radiative corrections - the leading logarithmic 
radiative corrections. The second problem is poor knowledge of charmonia DAs, which are the key ingredient of any 
calculation done within LC. 

Recently, the leadin g tw ist DAs of the 5-wave charmonia mesons have become the objects of intensive study 
[3 [3l|, US HI; HH I35L l36l. [37l l38j . Knowledge about these DAs allowed one to build some models for the S- wave 
charmonia DAs, that can be used in practical calculations. In this paper LC will be applied to the study the processes 
e+e~ — > J/^r] c , J/^rj' c , ip'r] c , ip'ri' c . It will be shown that with the models of DAs proposed in papers [H, [3J, [H[ LC 
predictions are in agreement with the results obtained at Belle and BaBar experiments. 

This paper is organized as follows. In the next sections a brief description of LC is given. In section III the formula 
for the amplitude of the process e + e~ — > VP, where V and P are vector and pseudoscalar mesons, is derived. In 
section IV this formula is used to calculate the cross sections of the processes e + e~ — > J/^r) c , J/^>r]' c , tp'r] c , ip'r]' c . The 
results of the calculation are discussed in section V. Finally, in the last section the results of this paper are summarized. 



II. BRIEF DESCRIPTION OF LIGHT CONE FORMALISM. 



In this section a brief description of light cone formalism (LC) will be given. As an example, let us consider hard 
exclusive process with single meson production. And let us assume that this meson is a pseudoscalar meson P. The 
presence of high energy scale Eh, which is of order of the characteristic energy of the hard exclusive process, allows 
one to apply factorization theorem for the amplitude of the process T 

T = Y / Cn x (P\O n \0), (1) 

n 

where the coefficient C n describes partons production at small distances, the matrix element (P\O n \0) describes 
hadronization of the partons which takes place at large distances. The sum is taken over all possible operators O n . 

For instance, the operators Q7 M 75<2, Qj^^D^D^Q, Qa^^GapQ are few examples of the operator O n . Actually, 
there are infinite number of the operators O n that contribute to the pseudoscalar meson production. 

The cross section of hard exclusive process can be expanded in inverse powers of the high energy scale Eh 

(TO (Tl , s 

CT = ™ + -^+1 + - ( 2 ) 

To determine if some operator contributes to a given term in 1/ Eh expansion one uses the concept of the twist of this 

operator (39|. Thus only the leading twist - the twist-2 operators Qzj$Q, Qz^^(zD)Q, Qzjc,(zD) 2 Q, contribute to 
the leading term in expansion l|2j). From this one sees that already at the leading order approximation infinite number 
of operators contribute to the cross section. Nevertheless, it is possible to cope with infinite number of contributions 
if one parameterizes all the twist-2 operators by the moments of some function 4>{x) 

(P(q)\Qz l5 (-iD tll )...(-W lln )Q\0) x z^...z^ = if P {qz) n+1 [ dx<j>(x)(2x - 1)", (3) 

Jo 

where q is the momentum of the pseudoscalar meson P, fp is the constant which is defined as (i 3 (?)|Q7^75<9|0) = 
ifpQ/jf: x is t ne fraction of momentum of the whole meson P carried by quark. The function (f>(x) is called the leading 
twist distribution amplitude (DA) . One can think of it as about the Fourier transform of the wave function of the 
meson P with lightlike distance between quarks. Using the definition of this function, factorization theorem ((TJ) can 
be rewritten as 

T= [ dxH{x)4>{x), (4) 
Jo 

where H{x) is the hard part of the amplitude, which describes small distance effects. This part of the amplitude 
can be calculated within perturbative QCD. As it was noted, the leading twist DA parameterizes infinite set of the 
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twist- 2 operators. This part of the amplitude describes hadronization of quark- antiquark pair at large distances and 
parameterizes nonpcrturbativc effects in the amplitude, ft should be noted that formula ^ resums the contributions 
of all the twist- 2 operators. If the meson P is a nonrelativistic meson, formula ([4} resums relativistic corrections to 
the amplitude T. 

Now let us consider radiative corrections to formula (|4]). The presence of two different energy scales, which are 
strongly separated Eh ~> Mp, gives rise to the appearance of large logarithm log E\jM P . This logarithm enhances 
the role of radiative corrections. The main contribution to amplitude ((J) comes from the leading logarithmic radiative 
corrections ~ (a s log E?/Mp) n . It turns out that these corrections can be taken into the account in formula (fj| as 
follows HUHJ 

T = / dxH(x, fi)(f)(x, [m). (5) 
Jo 

To resum the leading logarithmic corrections coming from all loops the scale \i should be taken of order of ~ Eh- 
The hard part of the amplitude H(x, /z) should be calculated at the tree level approximation. At this level H(x, //) 
depends on the renormalization scale \i only through the running of the strong coupling constant a s (/z). The rest of 
the leading logarithms are resummed in the DA <p(x,fj.) using renormalization group method (see Appendix B). It 
should be stressed that formula ((SJ) exactly resums the leading logarithmic radiative corrections. 

Commonly, to study the production of nonrelativistic mesons one uses effective theory NRQCD [25j . NRQCD deals 
with three energy scales toq 3> tuqv 3> ttiqv 2 , where tuq is the mass of the heavy quark Q, v <C 1 is relative velocity 
of quark antiquark pair. In the process of hard nonrelativistic meson production there appears one additional energy 
scale Eh which is much greater than all scales itiq , itiqv , itiqv 2 . Evidently, it is not possible to apply NRQCD at this 
scale. From the effective theory perspective, first, this large energy scale must be integrated out. And this is done 
through the taking into account renormalization group evolution of the DA <fi(x,fj,). 

This paper is devoted to the study of the process e + e~ — * VP, where V and P are vector and pseudoscalar mesons. 
Essential feature of this process is that for it the leading order contribution in I/-E7, expansion is zero [l4| • So, we deal 
with the next-to-leading twist process. For this process formula Q remains valid. The only difference is that now we 
have contributions coming from different twist- 2 and twist-3 DAs (see next section). As in the case of the leading twist 
process, formula Q resums relativistic corrections to the amplitude. However, if we consider the leading logarithmic 
radiative corrections to the amplitude, formula (O is incorrect. One can expect that some leading logarithms are lost. 
It is only possible to state that formula ([5]) resums the leading logarithms which appear in the amplitude due to the 
running of the a s and DAs. Below this approximation will be used. 

There is one common feature of all next-to-leading twist process. It is connected to the following fact: along with 
two particles twist-3 operators (for instance Q75Q) there appears operators of the type Q7 '/x7sG upQ ■ Evidently, these 
operators describe higher fock state \QQg > contribution to the amplitude of the process. NRQCD predicts that for 
nonrelativistic mesons such states are suppressed by higher powers of relative velocity of quark-aniquark pair inside 
the meson (25j . For this reason, in this paper the contribution of such states will be disregarded. 



III. THE AMPLITUDE OF THE PROCESS: e+e -> VP. 



In this section the amplitude of the process e + e~ — > VP, where V = J , ip' and P = rj c , r]' c will be considered. 
Two diagrams that give contribution to the amplitude of this process are shown in Fig [T] The other two can be 
obtained from the depicted ones by the charge conjugation. The amplitude of the process involved can be written in 
the following form: 

M = -4na HklMk2 \ v( Pl ,X)P(p 2 )\J-\0), (6) 
s p 

where a is the electromagnetic coupling constant, u(fci), life) are the electron and positron bispinors, s/s is the 
invariant mass of e + e~ system, J^ m is the electromagnetic current. The matrix element (V(pi, X)P(j>2)\J^ m \0) can 
be parameterized by the only formfactor F(s): 

(V( Pl ,X)P(p 2 )\j; m \0) = iq c F(s) e^M, (7) 

where q c is the charge of c quark, e\ is the polarization vector of the meson V{p\, A). The cross section of the process 
under consideration can be written as follows 

^e-^ VP) = ^(^)V(*)| 2 . (8) 
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FIG. 1: The diagrams that contribute to the process e + e — » V(pi, X)P(p2) at the leading order approximation in the strong 
coupling constant. 



In the last formula p is the momentum of the meson V in the center mass frame of the final mesons. 

To calculate the formfactor F(s) LC will be applied. As it was noted in the previous section, within LC the 
formfactor is a series in inverse powers of the characteristic energy of the process yfs. The leading order contribution 
to the formfactor F(s) in the 1/s expansion was obtained in paper [l4| . In derivation of this expression the authors 
of paper [l4[ disregarded the mass difference of the final mesons. For the mesons with different masses the expression 
for the formfactor F(s), which was derived in paper [l6j . can be written as follows 



\F(s) 



32tt 
~9~ 



f v fpMpM v 
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(9) 



dx 
o Jo 



1 MpZ t Z p V T (x,fi)Pp(y,Li) 
dy a s (n) J 



1 M c Z m Z t V T {x,n)PA(y,lj) 



M v 
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d(x,y) 
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1 

Mp 



d(x,y)s(x) Mp My 2 

1 (l-2y)V ± (x, f i)P A (y, f i) 
2M P s(y) d(x,y) 

(l + y)V A (x,n)P A (y,n) } 
d?{x,y) | ■ 



d(x,y) s(x) 



(10) 



(11) 



Where g 2 ~ (s 



M?, - Mp), x and y are the fractions of momenta carried by quark in the meson V and by quark 



in the meson P correspondingly, Pa, Pp, Vr, Vl, V±, Va are the DAs defined in Appendix B, My, Mp are the masses 
of the vector and pseudoscalar mesons correspondingly, M c — M^ is (ii — M* /s ), Z t and Z p are the renormalization 
factors of the local tensor and pseudoscalar currents, the dimensionless propagators d(x, y), s(x), s(y) are defined as 
follows: 



fc 2 

d{x,y) = — = 
% 
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a s (M c ) 



Zi 



a s (M 2 c ) 



Z„ 



2/(1 - y) 



» *(») = -3 



5 = 



x(l 



,(12) 




where cf = 4/3, b — 25/3, the definitions of the fourvectors fc, k x and k y are shown in Fig. [TJ 

It should be noted here that if one ignores renormalization group running, formula (J9j) is the analog of formula (|4|). 

So, in this case, it takes into account infinite series of the relativistic corrections to double charmonium production. 

If one takes into account the leading logarithmic corrections due to the DAs and the running of the strong coupling 

constant, the scale fi should be taken of order of the characteristic energy of the process ~ y/s. 

Before we proceed with the calculation of double charmonium production, the expression for the formfactor F(s) 

will be modified as follows: 
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1. Formula ([9]) was obtained at the first nonvanishing approximation in the 1/s expansion. This implies that 
all parameters in this formula must be taken at the same level of accuracy. For instance, in the expression ~ 
(s — My — Mp) last two terms are beyond the accuracy of calculation and these terms must be omitted. Below the 
following approximation will be used q$ — s. 

2. Now let us consider expressions (fT2j) for the propagators d(x,y), s(x), s(y). It is seen that all these expressions 
contain the terms proportional to ~ 5. These terms are of NLO approximation, so, according to the previous item, 
they must be omitted. Thus the propagators can be written as 

d(x,y)=xy, s(x) = x, s(y)=y. (13) 

However, if one substitutes these expressions to (TTU)) and takes into account the renormalization group evolution of 
the DAs, the divergence at the end point region (x, y ~ 0) appears. It should be noted that the motion of quark- 
antiquark pair in the end point regions x, y ~ 0, 1 is relativistic. So, if we consider the production of a nonrelativistic 
meson without the evolution of the DAs, the end point regions in the DA are strongly suppressed. For this reason 
the expression for the formfactor F(s) is free from the divergence. However, if the evolution of the DAs is taken into 
account, the divergence in the F(s) appears. In this paper this problem will be solved as follows. According to the 
definition of x 

EC+Pl (14) 



E M + pf ' 

where E c ,p^ are the energy and z component of the momenum of c-quark, E M ,p%* are the energy and z component 
of the momenum of meson M. If the energy and momenum of c-quark is much greater than it's mass, one can use the 
propagators in form (|13[) . since corrections to this approximation are suppressed. Let us now consider the kinematic 
region where c-quark is approximately at rest and x can be estimated as 

M* M* 



E M +pf 



M* here is the pole mass of c quark. In this region the accuracy of propagators (|13[) is not sufficient and one must 
take into account the corrections. It is clear that in any approach the corrections to the propagators regularize the 
whole expression for the formfactor F(s). This fact can be seen as follows: for any kinematical region of double quark 
and double antiquark production shown in Fig. 1 the squares of momenta of the gluon and quark propagators cannot 
be smaller than (2M*) 2 and (M* + Mv,p) 2 correspondingly. So, there is no divergence in the exact expressions for 
the propagators. 

In this paper this effect will be taken into account as follows: the propagators will be taken in form (|T3|) . but the 
integration in the expression for the formfactor F(s) will be done in the region x,y € {x m in, 1)- In other words to 
get rid of the singularity the cut off parameters x m i n is introduced. 

In principle, the calculation with propagators (TT2"]) is also possible. As it was noted already, in this case the terms 
proportional to the 8 play role of the regulator of expression (flQ)) . However, one can expect that the calculation done 
in this manner is less accurate and the result must be smaller than it is. To understand this one can apply the idea 
of duality of NRQCD and LC descriptions of the hard exclusive nonrelativistic mesons production: if LC expression 
for the formfactor F(s) is expanded in relative velocities of quark antiquark pairs of the mesons V and P , one will 
get NRQCD result for the formfactor. 2 

At the leading order approximation of NRQCD this idea can be reproduced if one takes infinitely narrow approx- 
imation for the DAs (Pi(x) = Vj(x) ~ S(x — 1/2), i = A, P, j = T,L,±,A ) and uses the following values of the 
parameters in formula (HJ) M v = M P = 2M*,M C = M*,fj, = M*,f VPT = {0\)s/M*. If this procedure is applied 
to (J9J) with propagators (|13|) . we will exactly reproduce the expression for the formfactor F(s) obtained within the 
leading order of NRQCD [J 

Fnrqcd = g g2 " (Oi)s, (16) 

where (Oi)s is the NRQCD matrix element Note that at the leading order approximation of NRQCD the 
formfactor F(s) scales exactly as 1/s 2 . So, within NRQCD 1/s 3 terms appear due to the relativistic corrections, 



2 It should be noted that these is no strict proof of this statement. However, one can expect that this statement is indeed true since the 
amplitude in NRQCD and LC can be expended in series of equivalent operators. Assuming that both theories can describe experiment 
one can expect that these expansions in both theories coincide. 
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which are suppressed as v 2 (v 2 is the characteristic velocity in charmonium) . In LC 1/s 3 terms can appear only due 
to the power corrections of the leading order result. So, applying the idea of duality of NRQCD and LC one can state 
that the power 1/s corrections to formula ([9]) with propagators (fl3|) are of order of ~ v 2 (M*) 2 / s . If we further apply 
the same procedure but with propagators (fT2"| . the formfactor F(s) will be different from the leading order NRQCD 
prediction 

F(s) = F NRQCD {1 - US + 0{5 2 )) = F NRQCD (l - +°(^)) • ( 17 ) 

To get the agreement with NRQCD prediction for the formfactor, one must expect rather large power correction to 
this result (14<5 x Fmrqcd) in LC so that to cancel the second term in the last equation. This correction appears at 
next-to-leading order approximation in 1/s expansion. So, the leading order approximation of LC with propagators 
([12)1 underestimates real result. It is not difficult to estimate the size of this effect using formula (fT7|) . For y/s = 10.6 
GeV and m* = 1.4 GeV the cross section calculated with propagators (Tl2"|) is smaller than that with propagators (fT3"|) 
by ~ 50%. Numerical calculation confirms this estimation. 

3. To calculate the cross section of the processes e + e~ — * J/^rj c , J /^r)' c , tp'r) c , tp'r)' c one needs to know the following 
constants: 

(J/y{p,e)\C la C\0) = fviMj/^ a (i/>'(p,e)\Cj a C\0) = iV 2 Afye Q , (18) 

(J/*(p,e)|Ca Q/3 C|0) M = ifri^ip ) (V ,/ (P ) e)|C f O'a/9C r |0>M = ifT2(n)(p a ep - Ppta), 

(Vc(p, e-)\C-f a j 5 C\0) = ifpip a (r]' c (p,e)\Cj ol 'y 5 C\0) = ifp2Pa- 

It should be noted that the operator Ca a pC is not renormalization group invariant. For this reason the constants 
/ti and /t2 depend on scale as 

c F 

fTM= {^)Y M ^ * =i ' 2 - (i9) 

In the derivation of formula ^ it was assumed that the tensor and vector constants fvi, fn(M c ) are connected to 
each other as fri{M c )/ fvi — 2M C /Myi, i = 1,2 [3, Gjj]. At the leading order approximation in relative velocity 
and strong coupling constant these relations are correct. However, they are violated due to radiative and relativistic 
corrections especially for the excited mesons. For this reason, below /y, fx will be treated as independent constants. 
Introducing the modifications described above, the expression for the formfactor F(s) can be rewritten as follows 



\F(s)\ = f 



f v fpMpM v 



Io , (20) 



l pi 



!),, 



/„=/ dxl dy aM { Ji_Mp_Z P VT(x,^)PpM _ h_ Mc_ Z m V T (x^)P A (y^) 
1 1 M v 2M C x 2 y M P 2M V x 2 y 



1 V L (x^)P A (y,fi) , (l-2y)Vj_(x,n)P A (v,li) , 1/ - 2M C \ 1 (l+y)V A (x,^P A (y,/ J L) 



+ -[l-f t Z n 



2M P x y 2M P xy 2 8 V M v ) Mp x 2 y 2 

where f t = /t(aO/ fv- Now we are ready to proceed with the calculation. 

IV. NUMERICAL RESULTS. 

To calculate the cross sections of the processes e + e~ — > J/^ij c , J/^ij' c , tp'r] c , tp'r]' c the following values of input 
parameters will be used: 

1. The strong coupling constant a s (/i) will be taken at the one loop approximation 

47T 

a «(^) = 75~i ( 2/A2V ( 22 ) 

Polog(M /A^J 

with A = 0.2 GeV, /3 = 25 /3. 

2. For the MS mass and the pole mass of c-quark the values M c = 1.2 GeV and M* = 1.4 GeV will be used. 
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H1H2 


&Exp X BrH 2 ^charged>2(R>) 


olo nrqcd{Po) 


o~nrqcd{Po) 


&poten model (f^) 


<TLc(fb) 




Bellof5] 


BaBar[6] 


m 


[2] 


[26] 




[21] 






25.6 ± 2.8 ± 3.4 


17.6±2.8±^ 


3.78 ±1.26 


5.5 


20.4 


17.61& 7 


22.2 ± 1.1 


14.41^ 


iP(2S)r] c (lS) 


16.3 ±4.6 ±3.9 




1.57 ±0.52 


3.7 






15.3 ±0.8 


10.41 9 7 ; 2 8 


V(lS)l7e(2S) 


16.5 ± 3.0 ± 2.4 


16.4±3.7±^ 


1.57 ±0.52 


3.7 






16.4 ± 0.8 


13.0lg;g 


iP(2S) Vc (2S) 


16.0 ±5.1 ±3.8 




0.65 ±0.22 


2.5 






9.6 ±0.5 





TABLE I: The second and third columns contain experimental results measured at Belle and Babar experiments. The 
BrH 2 ^charg<id>2 means the branching ratio of the decay of the hadron H2 into two charged particles. In the fourth and fifth 
columns the results of the leading order NRQCD obtained in papers P, 0] are shown. The NRQCD results obtained with 
inclusion of radiative and relativistic corrections [2r3. [2t1| are shown in columns six and seven. Potential model predictions [2ll ] 
for the cross sections are presented in column eight. Last column contains the values of the cross sections obtained in this 
paper. 

3. As it was noted in the previous section, to calculate the cross sections one needs constants (fT8|) . The constants 
fvi and fv2 can be determined directly from the experiment. The constants fxi and fx2 were calculated within 
NRQCD in paper [4(J. NRQCD formalism can also be used to determine the values of the last two constants fp\ and 
fp2 (see Appendix A). So, the calculation will be done with the following values of these constants 

0.173 ± 0.004 GeV 2 , /£ 2 = 0.092 ± 0.002 GeV 2 , (23) 

0.144 ± 0.016 GeV 2 , /^ 2 (M J/lt ) = 0.068 ± 0.022 GeV 2 , 
0.139 ± 0.048 GeV 2 , fp 2 = 0.068 ± 0.040 GeV 2 . 

4. The last inputs to formula ([20)) are the DAs. The models of the leading twist DAs will be taken from papers 
[33l [341. [3^. To build the models for the twist-3 DAs one can apply equations of motion. This procedure is described 
in detail in Appendix B. 

There are different sources of uncertainty to the results obtained in this paper. The most important uncertainties 
can be divided into the following groups: 

1. The uncertainty in the models of the distribution amplitudes 0i(cc,/x), which can be modeled by the variation of 
the parameters of these models (|B3[) . The calculation shows that for the processes e + e~ — > J/^r/ c , J/^>r)' c , ifi'Vc, ip'Vc- 
these uncertainties are not greater than ~ 5%, 12%, 28%, 40% correspondingly. It is seen that these uncertainties 
are not very large. This fact results from the property found in papers [331 ] : the evolution improves any model of DA. 

2. The uncertainty due to the radiative corrections. In the approach applied in this paper the leading logarithmic 
radiative corrections due to the evolution of the DAs and strong coupling constant were resummed. As it was noted 
above for the leading twist processes the resummation of the leading logarithms in the DAs and strong coupling 
constant is equivalent to the resummation of the leading logarithms in the whole amplitude. Unfortunately, this is 
not valid for the next-to-leading twist processes, for which some of the leading logarithms are lost. For this reason 
the radiative corrections should be estimated as ~ a s (^/s/2) log(s/4/(M*) 2 ) ~ 50%. 

3. The uncertainty due to the power corrections. This uncertainty is determined by the next-to-leading order 
contribution in the 1/s expansion. As it was noted above due to the application of the propagators in form (|13p one 
can hope that these corrections are suppressed by the square of relative velocity of quark antiquark pair in the meson. 
In the calculation these corrections will be estimated as ~ 4v 2 y M 2 ,/s ~ 25%. 

4. The uncertainty due to the regularization procedure. As it was noted above, to get rid of the divergence in 
the formfactor F(s) the cut of parameter was introduced. Evidently, our results depend on the value of the cut 
of parameter x m m • y/s, which is of order of the mass of c-quark (see formula (|15|) ). To estimate this source of 
uncertainty the cut off parameter is varied in the region x m i n ■ y/s = 1.0 — 1.6 GeV. The calculation shows that at 
x m in ■ \/s = 1.0 GeV the cross sections are increased by ~ 50% and at x m i n ■ y/s = 1.6 GeV the cross sections are 
decreased by ~ 20%. It should be noted that this source of uncertainty is closely connected with the uncertainty due 
to the radiative corrections. However, to understand this in detail one needs the theory which takes into account all 
the leading logarithmic corrections. 

5. The uncertainty in the values of constants \23\). It should be noted that this source of uncertainty is very impor- 
tant especially for the production of the excited states ip' , rj' c . Thus, for the processes e + e~ — * J/^n c , ip'vc, J/^n^, ^}'v( c 
the errors due to the uncertainties in the values of constants ([23)) are ~ 35%, 40%, 60%, 65% correspondingly. 

Adding all the uncertainties in quadrature one gets the total errors of the calculations. 



fh 
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H\H2 


OLO NRQCD 




otot (fb) 


4>(lS) Vc (lS) 


1.9 


6.3 


14.4 


i>(2S)r) c (lS) 


1.0 


6.2 


13.0 


V(lS)r? c (25) 


1.0 


7.8 


10.4 


i>(2S)r) c (2S) 


0.53 


7.2 


9.0 



TABLE II: The second column contains the values of the cross sections obtained at the leading NRQCD approximation. The 
third column contains the cross sections calculated in the following approximation: all relativistic corrections are resummed, 
but the leading logarithmic radiative corrections are not taken into account. The last column represents the results obtained 
if the relativistic and leading logarithmic radiative corrections to the amplitude are taken into account simultaneously. 

The results of the calculation are presented in Table I. The second and third columns contain experimental results 
measured at Belle and Babar experiments. In the fourth and fifth columns the results of the leading order of NRQCD 
approach obtained in papers [3, 0] & r e shown. The NRQCD results obtained with inclusion of radiative and relativistic 
corrections 0, [27j are shown in columns six and seven. Potential model predictions [2l[ for the cross sections are 
presented in column eight. Last column contains the values of the cross sections obtained in this paper. 

V. DISCUSSION. 

It is seen from Tab. I that within the accuracy of the calculation the results of this paper are in agreement with 
Belle and BaBar experiments. It is also seen that the uncertainty of the calculation is rather large. There are two 
very important sources of uncertainty. The first one is the theoretical problem with taking into account of all leading 
logarithmic radiative corrections to the amplitude of the next-to-leading twist processes. It can be estimated as 
~ 50 — 70% of the cross sections. This source of uncertainty can be reduced if the theory, which takes into account 
all leading logarithmic corrections, is created. The second very important source of uncertainty is poor knowledge of 
constants (IT51) . For some processes this uncertainty can reach 60%. The values of these constants used in this paper 
can be considered as the first estimation of their real values. So, theoretical and experimental study of these constants 
can greatly improve the accuracy of any predictions done within LC. 

Next, let us discuss the results obtained in other papers and compare them with the results of this paper. First, let 
us consider the results of the leading order NRQCD predictions [HQ shown in columns four and five of Tab. I. These 
results are approximately by an order of magnitude smaller than the cross sections measured at the experiments. At 
the same time LC predictions are in reasonable agreement with the experiments. This facts lead to the question: 
why LC predictions are much greater than the leading order NRQCD predictions? The answer to this question can 
be given within LC. As it was noted in section III, LC can reproduce the leading order NRQCD results. To do this 
the renormalization group evolution of the constants and DAs will be disregarded and all parameters will be taken 
at the leading NRQCD approximation: the constants fxi, fpi are equal to the constant ,fv%,i = 1,2, which can be 
determined from the leptonic decay width, the mass M c = M* = 1.4 GeV, My = Mp = 2M*, all DAs are taken at 
the infinitely narrow approximation ~ 5(x — 1/2). Thus one gets the results shown in the second column of Tab. II. At 
the second step all parameters used in the calculation are taken at their central values, but without renormalization 
group evolution. At this step infinite series of the relativistic corrections are resummed, but the leading logarithmic 
radiative corrections are not taken into account (see section II). The results are shown in the third column of Tab. 
II. At the last step, renormalization group evolution is taken into account and the results are presented in the last 
column of Tab. II. Within LC this means that the relativistic and leading logarithmic radiative corrections are taken 
into account simultaneously. 

From Tab. II one sees that the relativistic and leading logarithmic radiative corrections taken into account si- 
multaneously dramatically enhance the leading NRQCD predictions and bring the agreement with Belle and BaBar 
experiments. Very important conclusion which can be drawn from this result is that in hard exclusive processes rela- 
tivistic and leading logarithmic radiative corrections play very important role and the consideration of such processes 
at the leading NRQCD approximation is unreliable. 

Looking to the results of Tab. II one can also draw a conclusion that relativistic corrections alone cannot describe 
the experimental results. This conclusion is in agreement with the results of papers [j], [ID, [H, H(| ■ In these papers the 
cross section of the process e + e~ — * J /^r] c is enhanced due to the relativistic corrections by a factor of ~ 2 — 3 what 
is in agreement with the results obtained in this paper. At the same time this conclusion is in disagreement with the 
results of paper [2l[ (see Tab. I column eight), where the authors tried to attribute the disagreement between theory 
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and experiment only to the relativistic corrections, which were calculated within potential model. From the results 
shown in Tab. II one sees that this approximation is realistic only for the production of excited states e + e~ — > ip'v'c- 
In this case the relativistic corrections are much more important than the leading logarithmic radiative corrections. 
The authors of papers 

QUI 

took into account only the part of the leading logarithmic radiative corrections which 
appears due to the evolution of the constant fx and the running mass of c-quark. The evolution of the DAs was 
disregarded. The calculation done in this paper shows that this approximation can be applied only to the DAs of 
the 2S charmonia mesons, for which the evolution is not very important ( see paper [35j). As it was shown in papers 
[33l l34j the evolution of the IS state charmonia DAs is very important. To compensate the effect of the evolution of 
the DAs the authors of paper [3] proposed rather wide model of the DAs with so called relativistic tail 3 (see paper 
|3ll]). For instance, the characteristic velocity of the IS* charmonia (see formula (|A2[) ) with this DA can be estimated 
as 



(Ois ~ 3(Ois = 0-39, (24) 



which is much larger than (v 2 )is — 0.2 — 0.3 calculated in papers [191. 133, l34l. l4lT |42|| . 

At the end of this section let us consider how the disagreement between theory and experiment can be resolved 
within NRQCD. The authors of papers 0, [2tJ took into account the relativistic and one loop radiative corrections 
and got the following values of the cross section of the process e + e~ — > J/^ijc 

a(e + e~ -► J/^r/c) = 20.4 fb HI, 

cr(e + e~ -► J/^r/c) = 17.6t8° 3 7 ft> S3- (25) 

So, similarly to this paper, the authors of papers [26|, [2tJ resolved the disagreement through the taking into account 
of the relativistic and radiative corrections. 

Central values ([23)) are larger than the central values of the cross section obtained in this paper. Below it will be 
shown that central values (|25j) are overestimated. A simple way to understand why the cross section was overestimated 
is to consider paper which has similar problem. In this paper explicit expression for the one loop radiative 
corrections to the cross section of the process e + e~ — > J /^rj c were calculated. The result of this paper can be written 
in the form 

a = a"x^ + ^M K y (26) 

where cr° is the cross section at the leading order approximation, the expression for the factor K can be found in [ill ]. 
The cross section a is proportional to \Rj/^,(0)\ 2 x \R Vc (0)\ 2 , where Rj/&(r), Rq c {r) are the radial wave functions 
of the J/ty and r\ c mesons. It should be noted here that the cross section a is very sensitive to the values of the 
wave functions at the origin Rj/^(0), R r)c (0), so it is very important how these parameters were calculated. In the 
calculation the authors of [ll[ took |i?j/^,(0)| = |-R?7 C (0)| and the value of the |i?j/^(0)| was taken from the leptonic 
decay width T ee 

l«,/*(»>l 2 - r^^T W 

3 7T ^ c 

Now few comments are in order. First, this expression is taken at the next-to-leading order approximation in the 
strong coupling constant, what means that some part of the one loop radiative corrections is put to the a . Second, 
if formula (|27p is expanded in the a s , one will get infinite series in the strong coupling constant. So, the application 
of formula (I27j) in the calculation of the wave function at the origin is equivalent to the statement that one knows all 
infinite series of the radiative corrections to the wave function. Evidently, this is not correct. To be consistent with 
the one loop approximation applied in paper [ll[ , one should expand expression (|27|) in the strong coupling constant 
and than leave only the first term in this expansion 

i«4 + f^)%^. m 



3 It should be noted that relativistic tail of the DAs used in this paper is absent at the scale fi ~ M* and appears due to the evolution of 
the DAs at larger scales. 
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It is not difficult to see that the \Rj/^(0)\ 2 calculated from formula ((27)) is greater than that calculated form formula 
(|28p by a factor 1/(1 — (16a s /37r) 2 ). The cross section calculated using formula (|2"7| is greater than the cross section 
calculated using formula ([28]) by a factor 1/(1— (16a s /37r) 2 ) 2 ~ 1.5 for a s = 0.25. So, the values of the wave functions 
at the origin calculated in paper [ll| were overestimated, what led to the overestimation of the cross section by a 
factor of ~ 1.5. 

Similar overestimation of the wave functions at the origin and, as the result, overestimation of the cross section 
takes place in papers [26|, [27[. For instance, in paper (2?} the authors used the value of the NRQCD matrix element 
(0i)j/*, which is analog of the wave function at the origin, calculated in paper [42| using the formula 

1 3M 2 /q ,T ee 
(l-/((f 2 ),;/*)-|v) 2 8na 2 qi 

where (v 2 )j/* is defined in (|A2[) . f(x) = xj (3(1 + x + y/1 + x)). Similar expression can be written for the (Oi) ))c ■ It 
is clear that the application of formula (|29p is equivalent to the statement that one knows not only full a s corrections 
in front of the leading NRQCD operator (Oi)j/* but also full a s corrections in front of all operators that control 
relativistic corrections (i> n )j/$ to the amplitude of the leptonic decay. Evidently, this is not correct. To make formula 
(|29f more consistent with the approach applied in paper [27}, it should be modified as follows 

7-2 



SMj/yTee / 1 16 a t 



<°i> J„ = ( r, 777^ ~* + ~ ) ■ (30) 



In last formula all problems mentioned above are resolved. Similarly, one can improve the expression for the operator 
{0\)n a - With numerical values of paper [12] the factor of the overestimation is ~ 1.5. Taking into account this factor, 
the cross section 17.6 fb obtained in paper [13] should be changed to 11.7 fb. 

Similar analysis can be done for the result of paper [26j . The calculation shows that the cross section is overestimated 
by a factor of ~ 1.5. Taking into account this factor the central value of the cross sections now is ~ 13.6 fb. 



VI. CONCLUSION. 



This paper is devoted to the study of the processes e + e~ — > 3 j^r\ Cl J /^>i]' c , ip'rj c , ip'v'c within Light Cone Formalism 
(LC). The amplitude for these processes was first derived in paper [14j |. In the present paper the amplitude was 
modified so that to archive better accuracy of the calculation and to resolve some questions raised in [3l|. 

To calculate the cross sections of the processes under study one needs the twist-2 and twist-3 distribution amplitudes 
(DA). The models of the twist-2 DAs of the S- wave charmonia were proposed in papers [H, [H, [H[ . To get the twist-3 
DAs, equations of motion were applied. It turns out that if one ignores the contribution arising from higher fock 
states, the twist-3 DAs can be unambiguously determined. 

Using the models of the DAs, the cross sections of the processes e + e~ — > J/^r} c , J/^T]' C , ip'r) c , vp'v'c were calculated. 
Within the error of the calculation the results of this study are in agreement with Belle and BaBar experiments. 
In addition, the question - why LC predictions are much greater than the leading order NRQCD predictions - was 
studied. Numerical results of the calculation shows that large disagreement between LC and the leading NRQCD 
predictions can be attributed to large contribution of relativistic and radiative corrections. From this results one can 
draw a conclusion that in hard exclusive processes relativistic and radiative corrections play very important role and 
the consideration of such processes at the leading NRQCD approximation is unreliable. 

The results of this paper are in agreement with recent NRQCD study of the process e + e _ — > J j^r\ c [H, H3| where 
the authors took into account relativistic and one loop radiative corrections. However, in the present paper it was 
shown that the results of these papers are overestimated by a factor 1.5. 
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APPENDIX A: CALCULATION OF THE CONSTANTS /pi,/p 2 . 

To calculate the values of the constants fn, i = 1, 2 one can apply NRQCD formalism. At the NLO approximation 
of NRQCD the constants fpi can be written as follows [25|, l43j 



where 



flM? = ^k + .\|0)<0|x + y -|i>,} (l-4^-(A)- (AD 



1 (0|xW^) (A2) 
{V)l ~ ml <0| X +^> ' (A2) 

Pi here is the ?y c meson if i = 1 and the 77^, meson if i = 2. To continue the calculation one also needs NRQCD 
expression for the decay width T [rj c — > 77] : 

, 1 , . / 7r 2 — 20 a s (v 2 )„ \ , . . 

T[r lc ^ 11 } = ^-( Vc \^+ x \0)(0\ x + ip \r ]c ) x h +___!- l^j (A3) 

To determine the constant fp\ let us express this constant through the width T [rj c — > 77] 4 

where M* is the mole mass of c-quark. The value of the constant fpi will be calculated with the following set of 
parameters: V [ry c -> 77] = 7.2 ±0.7 ±2.0 [44], a s = 0.25, (v 2 ) Vc = 0.25 41], (u 2 )^/ = 0.54 [35], M* = 1.4 ±0.2 GeV. 
To estimate the error of the calculation one should take into account that within NRQCD the constant is double series 
in the relativistic and radiative corrections. At the NNLO approximation one has the relativistic corrections ~ (v 2 ) 2 , 
the radiative corrections to the short distance coefficient of the operator {Se)ip\Vi(e)) ~ a 2 and the radiative 

corrections to the short distance coefficient of the operator (0\x + {<Je) CD) 2 ip\Vi(e)} that can be estimated as ~ a s (v 2 ). 
In addition, there is an experimental uncertainty in the measurement of T [rj c — » 77] and the uncertainty in the m c . 
Adding all these uncertainties in quadrature one can estimate the error of the calculation. Thus one gets 

fli = 0.139 ± 0.048 GeV 2 . (A5) 

Unfortunately, one cannot apply formula (IA4[) to get the value of the fp2- Since today only the product T [rf c — > 77] x 
Br(rj' c — * KgK^TT ) has been measured [44J and there is no model independent way to determine T [r]' c — > 77]. 

To estimate the value the fp2 one could use the value of the constants fv2 of /t2(M//$), which are equal up to 
the relativistic correction and radiative corrections. In this paper the constant /t2(AO/*) Wl U be taken to get the 
value of fp2 

fp 2 = 0.068 ± 0.040 GeV 2 . (A6) 

APPENDIX B: MODELS FOR THE DISTRIBUTION AMPLITUDES. 

In this section models and renormalization group evolution of the DAs needed in the calculation will be considered. 
The DAs of the vector meson V are defined as follows [l4| 

(V(p, e)\Q a (z)Qp(-z)\0)n = f -^- J dxjW**- 1 ) { (t - V±(x, n)+p^-V L (x, M ) + 

Mfl) 1 «^P» Vr(x,/x) + h ^ e w /^V V A ( X ,A . (Bl) 



fv M v 



In the derivation of J A4I ) the expression M Vc = 2M* + M*(v 2 ) Vc was used. 
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The DAs of the pseudoscalar meson are defined as 



(P(p)\Q a {z)Q fl (-z)\0) ll = i 



JpM v 



(B2) 



/3a 



where x is the fraction of momentum of the meson V carried by quark, the constants /t(m), fv are defined in equations 
(HHJ), M c (/x) is the running mass of c-quark. Expressions (IB 1 1 IB2 j) are defined at the scale fx. 

The models of the leading twist DAs Vl(x, fi),Vr{x, ^),Pa{x, fi) were proposed in papers [H, [H, [35| . According to 
these models the functions Pa{x,/j, ~ M c ),Vr(x, fi <~ M c ),Vl(x, fi ~ M c ) are equal to the function <j>\s{x,U ~ M c ) 
for the IS" state mesons and to the function <^>2s(x, A* ~ -^e) f° r the 25 state mesons, which have the form |30| 



fcs(x,V ' 



Mr 



lfias{x) Exp 

< y 9 as (x)(a + 



4x(l - x) 



2x - 11 



Exp 



4x(l - x) 



(B3) 



where ip as (x) = 6x(l — x) is the asymptotic function. For the IS charmonium states the constant b can vary within 
the interval 3.8 ± 0.7. For the 2S charmonium states the constants a and b can vary within the intervals 0.03+q'q?, 
and 2.5^Q g correspondingly. The renormalization group evolution of the leading twist DAs is well known [30l | and it 
can be written in the form 



(x, /it) = 6x(l — x) 



n=2,4. 



Mm) 

a s (Mo) 



,/bo 



a„(/io)C I 3 / 2 (2x- 1) 



(B4) 



where a n (fJ,o) is the coefficient of the expansion in Gegenbauer polynomials C^J 2 {z) at scale /x , the constant &o = 25/3, 
the anomalous dimensions e n for the functions Vl, Pa are defined as 



(ra + l)(ra + 2) 



n+l 



J=2 



(B5) 



the anomalous dimension e„ for the function Vp 



= -U 



n+l ^ 



J'=2 



(B6) 



Further let us consider the twist-3 DA Pp(x,fi). It turns out the if one ignores higher fock states it is possible to 
connect the twist-3 DA Pp(x,fi) to the twist-2 DA Pa(x, fi) using equations of motion [45[ 



(Dp = $nO + 



(B7) 



where £ = 2x — 1, (£™)p,a are the moments of the DAs Pp and Pa correspondingly, r(ji) = 4M c (/i) 2 /M 2 . To solve 
these equations one can expand the Pp in a series of Gegenbauer polynomials GI/ 2 (z) |4(| 

Pp(x,fi) 



1+ £ 6 n (/i)Cy2(2x-i; 

n=2,4.. 



(B8) 



Substituting expressions for the DAs Pp and (|B4[) and (|B8|) to equations of motion (|B7[) , one can solve these equa- 
tions recursively and relate the coefficients b n (fi) to the coefficients a n (fi). For instance, for the first tree coefficients 
one can get the following formulas 

5 27 81 13 39 195 

b 2 (n) = -~r{n), b 4 (p) = --r(n) - —r a 2 (n), b 6 (fi) = -— r(fi) - — r(/x)a 2 (/x) - — rO)a 4 (»' (B9) 

It is clear that relativistic motion in the nonrelativistic system must be strongly suppressed. The behavior of DA in 
the end point region x ~ 0, 1 is determined by the relativistic motion. So, the DA of nonrelativistic system must be 
strongly suppressed in the end point region. As it was shown in paper (33j for the leading twist DAs such suppression 
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can be achieved if there is fine tuning of the coefficients a n (fi) at the scale fi ~ M c . Similarly, to get the suppression 
in the end point region for the function Pp one requires the fine tuning in the coefficients b n {M c ) 1 which can be 

2 

achieved through the fine tuning of the constants r(M c ) and a n (M c ). If we put r(M c ) — AM C /M 2 the fine tuning 
for the higher moments will be broken what leads to large relativistic motion in the end point region. To avoid this 
problem it will be assumed that r(M c ) is a free parameter, which will be adjusted through the requirement that the 
moments of the Pp(x,fi ~ M c ) must be equal to the moments of the Pa(x,/j, ~ M c ) to the leading order in relative 
velocity expansion. Using the last statement and (|B7[) one can obtain the expression for the r(M c ) 

(BIO) 



r(M c ) = 1 - 3(e 2 >A + 9((C 2 )a) 2 - - 27((£ 2 }a) 3 + W^UifU - 7(f) a + 0(v 8 ). 



All moments in this expression are taken at scale fj, ~ M c . The exact expression for the r(M c ) can be written in the 
following form 



(m c ) y_i ^(i-e) 



;Pa 



(Bll) 



where £ = 2x — 1. 

The same approach can be applied to the twist-3 DAs Va,V±. These functions can be expanded in a series of 
Gegenbauer polynomials (471 ]: 



V±(x,fi) = 
Va{x, [i) = 6x(l — x) 



1+ E c n (^)C 1 J 2 (2x-l) 

n=2,4.. 

1+ E d n (v)C?J 2 (2x-l) 



(B12) 



The coefficients c n (/z) and d„(/i) can be related to the coefficients a n (p) of the functions Vl, Vt through the equations 
of motion [H 



(n + = <Dl + ^-il(l - 5(»))(C- 2 )a, 

i(n + 2)(1 - <5(m))(CU = (Di - <Km)(Ot, 



(B13) 



where (^") l,t,±,a are the moments of the DAs Vl,Vt,V±,Va, S(/j,) = 2fx{n) / fvM c {n) / My ■ Similarly, to the constant 
r(M c ) the value of the constant S(M C ) must be adjusted. We adjust it through the requirement that the moments of 
the DA Vj_ must be equal to the moments of the Vl to the leading order approximation in relative velocity expansion. 
Thus one can get 



S(M C ) = 1-1 2<£ 2 )l - 6((£ 2 ) L ) 2 + 4(C 4 )l - 22<e 2 ) L (£ 4 )L + 18«e 2 )L) 3 + 6<£ 6 )l + 0(v 8 ) 



(B14) 



It should be noted that with this value of the S(M C ) the moments of the function Va are equal to the moments of the 
Vl to the leading order in relative velocity. The exact expression for the 8{M C ) is 



6(M C ) 



(B15) 
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